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satellite  about  LA  is  per  io  ruled.  A  p  r  ■  ;>  o  s  l  two  J  iia  us  i  ona  1 

very  restricted  orbit  is  used  to  supply  the  initial 
conditions  required  for  the  search.  An  ophom.eMs  of  high 
accuracy  is  generated  from  a  specific  date  and  tine  using 
actual  positions  for  the  sun  and  moon.  The  generated  sun  and 
moon  position  and  velocity  vectors  are  used  in  the 
integration  of  the  system's  equations  of  motion.  A  stable 
orbit  is  found  and  is  tested  for  its  length  of  stability. 
The  orbit  is  found  to  have  a  stable  lifetime  in  excess  of  six 
hundred  lunar  synodic  months.  The  sensitivity  of  the  orbit 
to  the  sun's  and  moon's  position  is  tested  and  found  to  be 
only  slightly  sensitive  for  an  error  in  position  of  one 
quarter  day.  Finally,  a  predicted  180°  out  of  phase  orbit  is 

found  and  is  determined  to  be  only  marginally  stable. 
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and  Earth-based  laser  weapons. 

Studies  undertaken  so  far  have  largely  dealt  with  two 
dimensional  analysis  of  the  very  restricted  four  body 
problem.  Some  of  the  restictions  have  been  circular  orbits 
for  both  the  moon  and  sun,  an  unrealistic  mass  for  the  sun, 
and  planar  motion  for  all  bodies.  Some  analysis  has  been 
undertaken  in  the  past  to  study  the  three  dimensional 
problem,  but  these  have  always  been  restricted  in  their  scope 
and  lead  to  a  vague  conclusion  about  the  existence  of  three 
dimensional  orbital  stability  in  the  Lagrangian  vicinities. 
This  report  will  attempt  to  show  actual  stability  about  L4 
for  a  period  of  at  least  fifty  years,  ana  to  set  conditions 
for  further  studies  that  follow  this  report. 

Most  recently,  at  the  Air  Force  Institute  of 
Technology,  where  this  report  was  written,  two  studies  were 
completed,  one  by  Major  William  Beckman  (Kef  3),  and  the 
other  by  Captain  John  Wheeler  (Ref  7).  Each  did  a  study  of 
two  dimensions  of  orbital  stability  about  L4.  Both  also 
lacked  an  analysis  of  the  moon  and  sun's  actual  positions  and 
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proposed  orbit.  iiis  orbit  was  Linn  Leu  by  the  .'a  r  ions 
constraints  that  were  imposed  on  the  problem.  The  major 
limitation  came  about  by  assuming  the  moon  and  sun  to  be  in 
circular  orbits  a  b  o  ■.  t  their  respective  system  barycenters. 
Capt.  Wheeler's  conclusions  at  the  end  of  his  study 
indicated  linear  stability  exists  for  his  system. 

Major  Beekman's  study  was  based  on  three  reports,  one 
of  them  being  Capt.  Wheeler's.  The  other  reports  were  by 
Kolenkiewicz  and  Carpenter,  in  1968,  (Ref  6),  and  by  Barkham, 
Modi,  and  Soudack  in  1975,  (Ref  2).  In  his  investigation. 
Major  Beekman  confirmed  the  Wheeler  model's  stability,  for  a 
period  of  at  least  twenty  years,  by  removing  the  restriction 
of  circular  orbits.  he  also  showed  the  utiier  u  r  d  i  L  s  studied 
were  marginally  stable  in  the  same  time  frame,  even  though 
the  planar  restriction  was  still  in  force. 

Problem  and  Scope 

The  search  for  a  three  dimensional  stable  orbit  about 
L4  and  L5  is  a  problem  of  increasing  importance  in  the 
evolution  of  space  exploration.  For  any  serious  long  term 
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dimensionality  and  the  unperturbed  circular  orbits  of  the 
moon  and  sun  . 
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I  1  . 


In  l  h  e  analysis  of  l  h  o  probk-i.:,  the  equation  s  of 
motion  of  the  moon  and  the  colony  are  those  given  by  T.  A  . 
Heppenheimer,  in  an  article  published  in  197  8,  (Ref  5).  i  i  s 

equations  are  given  in  rectangular  nonrotating  Earth-centered 
inertial  coordinates  for  two  dimensions.  These  equations  are 
then  simply  expanded  to  three  dimensions.  The  sun's  motion 
is  then  described  in  the  same  coordinate  system  as  a  two  body 
problem.  Since  the  equations,  are  given  in  only  two 

dimensions  the  expansion  to  three  dimensions  follows  the 
derivation  of  the  two  dimensional  case. 

The  equations  of  motion  for  a  general 
problem,  in  rectangular  coordinates  is  given  by: 

n 

ht.i  (?  -F  ) 

J  i _ i 

3 

rij 

where  r  =|r  -r  |.  Letting  i=l  be  the  reference  body,  the 

i  j  1  j  i 

Earth,  i-2  the  body  whose  motion  we  wish  to  study,  and  i=3 
and  i**4  be  the  indices  of  the  two  other  bodies  in  the  system. 
For  the  Earth-Sun  system: 


"n-body" 


i;g  li-l 
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From  EQ  II-3,  the  equations  of  the  motion  of  the  moon  are 
written,  also  realizing  the  satellite  has  no  effect  on  the 
moon ' s  motion: 
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EQ  1 1 - 6 

Using  the  Earth  as  a  convenient  reference  frame,  the 
equations  finally  become: 


EQ  I  I - 9 

The  equations  of  motion  for  analysis  purposes  are 
then  represented  in  state  vector  form  for  ease  in  handling. 
See  Appendix  A  for  the  subroutine  pertaining  to  the  equations 
of  motion. 
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Ref  8.  Since  tlie  moon 
right  ascension,  a, 
must  be  transformed  to 
rectangular  coordinates 


racing,  rectangular  coordinate  system 
The  position  vectors  associated 
-  sun  in  this  frame  are  developed  from 


vector 

given 

by 

Ref 

8  is  described 

i  n 

dec  1  in 

ation, 

d  , 

and 

Earth  radii,  r 

.  i  t 

the 

frame 

we 

wish  to  use. 

The 

are  g 

i v e n  by 

x=r*cos(a)*cos(d) 

y*r*sin ( a) *cos ( d )  EQ  11-10 

z-r*sin(d) 

# 

and  the  velocity  vector  elements  are: 

x*‘r*cos(a)cos(d)-r*a*sin(a)cos(d)-r*<5*cos(a)sln(d) 

y«r*sin(a)cos(d)-r*a*cos(a)cos(d)-r*d*sin(a)sin(d) 

z«r*sin(d)+r*d*cos(d) 

EQ  1 1 - 1 1 

The  vectors  are  still  not  in  the  proper  frame  and  need  to  be 
rotated  to  the  ecliptic.  If  e  is  the  obliquity  of  the 
ecliptic*  then  the  transformation  matrix  for  this  is: 
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0 


i  n 

0  r  o  :'.(•■  '  s  i  n  <  •  ' 

0  -sin(e)  cos(e) 


EQ  11-12 


The  sun's  position  sector  is  already  in  Earth-centered 
rectangular  coordinates  and  only  needs  to  be  rotated  to  the 
proper  frame  by  the  use  of  the  above  transformation  matrix. 


The  frame  for  the  analysis  of  the  problem  will  be  an 
Earth-centered  ecliptic  nonrotating  rectangular  system.  The 
X-axis  will  point  toward  the  vernal  equinox  and  the  Z-axis 
will  be  perpendicular  to  the  ecliptic  having  the  XY-plane 
coincident  with  the  ecliptic  plane.  The  frame  for  the 
presentation  of  the  output  of  the  analysis  will  be  a  rotating 
frame  with  the  x-axis  through  the  center  of  mass  of  the  moon. 
It  also  will  be  an  Earth-centered  rectangular  ecliptic  frame. 
See  Fig  1  and  Fig  2  for  a  pictorial  representation  of  each 
coordinate  system. 

Ephemer is  Generation 

Ref  8  provides  the  position  vector  for  both  the  moon 
and  sun,  but  does  not  provide  a  velocity  vector  for  each 
appropriate  time  step.  In  order  to  integrate  the  equations 
of  motion,  the  velocities  must  be  known  at  any  given  time. 
Since  Ref  8  provides  position  at  a  known  time,  a  central 
difference  velocity  can  be  determined.  This  velocity  is 
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Using  the  position  and  velocity  vectors  as  initial 

conditions,  the  equations  of  notion  are  integrated  forward  to 

a  reference  time,  t  ,  where  the  position  vector,  p  ,  was 

o  o 

known.  The  reference  time  selected  was  a  function  of  the 
time  step  available  from  Ref  8,  ten  days  for  the  sun  and  one 
half  hour  for  the  moon.  The  corrector  was  first  applied  only 
to  the  sun's  velocity  and,  after  repetition  at  longer 
reference  periods,  the  corrector  was  applied  to  the  moon's 
velocity.  The  outcome  of  the  corrector,  the  velocity  vectors 
of  the  moon  and  aun  at  the  initial  time,  was  substituted  for 
the  crude  velocities  of  Ref  8. 

The  start  date  for  the  ephemeris  generation  was  chosen 
as  5  Jan  1979,  of  the  Equinox  of  1950.  This  particular  date 

was  chosen  because  the  moon  was  relatively  near  the 
equatorial  plane  of  the  Earth.  Using  a  time  when  the 
Z-component  of  the  moon's  velocity  vector  approaches  zero 
will  give  highly  inaccurate  velocities  when  computed  by 
central  differences  and  will  require  more  repetitions  of  the 
velocity  corrector  to  achieve  the  velocity  desired. 
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The  first  sUp  ni  c  o  m  put  i  :i ,  t  i  i  «.•  .  i  l.u'ity  corrector  is 
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o  1 

is  compared  to  tile  r  e  t  e  r  e  u  c  e  position,  p  ,  and  this  vector, 

o 

dr  =  p  -p  ,  is  stored  for  1 ..  i  e  r  n  -  pet  ioit.  of  mot  i  o 

L  o 

arc-  then  into  rat.  i  tor  .ar.f  tr.ree  ore  t  i..iea  from  ti .  r  t  :  >.  1 

conditions  with  the  initial  position  vectors  ano  tile 
following  velocity  vectors  in  turn: 


r  — 

r 

-  ■■ 

x+dv 

X 

X 

y 

» 

y+dv 

y 

z 

z 

z+dv 

EQ  11-13 


where  dv  is  a  small  velocity  increment.  It  should  be  noted 

that  two  of  the  velocity  components  are  the  original  initial 

velocity  conditions  for  each  integration.  Each  time  the 

equations  are  integrated  forward  a  new  columnar  position 

vector,  p  ,  is  formed.  Subtracting  p  from  p  ,  one  obtains 
i  o  i 

the  vector  dp  =p  -p  ,  where  i=l,2,3.  After  the  three 
i  i  o 

integrations  these  position  vectors  may  be  combined  into 

matrix  form,  A-[dp  :dp  :dp  ].  Dividing  this  matrix  by  the 

12  3 

delta  velocity  that  was  added  to  the  initial  velocity 
vectors,  a  differential  matrix  is  obtained  of  the  form: 
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.  This 

velocity  corrector 

vector  is 

th 

e  n  added 

t  0 

the 

initial 

conditions  velocit 

v  vector. 

This 

procedure 

is 

then 

repeated 

using  a  new  reference  time,  until  the  velocity  vector  yields 
a  postion  vector  at  the  end  of  the  ephemeris  span  to  within 
the  accuracy  desired.  See  Appendix  B  for  a  subroutine 
pertaining  to  these  calculations. 


Constants 


The  primary  constants  used  for  the  problem  analysis 
were  obtained  from  Ref  [8:529].  All  secondary  constants  were 
derived  from  these  values.  The  accuracy  of  the  constants  in 
Ref  8  is  on  the  order  of  six  digits  or  less  for  masses  and 
distances  but  this  will  be  shown  to  be  adequate  for  the 
problem  analysis.  The  use  of  constants  other  than  those 
obtained  in  Ref  8  are  used  only  in  the  duplication  of  the 
Wheeler  model  and  are  obtained  from  Ref's  3  and  7. 
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The  Wheeler  Model 

The  Wheeler  orbit  should  be  obtainable  from  the  method 
of  analysis.  In  order  to  achieve  this,  the  equations  of 
motion  are  modified  to  place  the  moon  and  sun  in  circular 
orbits  about  the  Earth.  Wheeler's  constants  are  then 
corrected  to  the  unit  constants  described  in  the  preceeding 
section.  The  initial  conditions  are  then  transformed  to  the 
frame  coordinates  being  used,  and  the  equations  of  motion  are 
then  integrated  forward  for  the  appropriate  time  span,  one 
lunar  synodic  month.  See  Fig  3  for  a  representation  of  the 
Wheeler  system  and  Fig  4  for  the  predicted  Wheeler  orbit. 

The  initial  conditions  .re  selectee  i  r  o  ..  the 

restrictions  of  the  Wheeler  model.  The  sun  and  moon  are  in 
circular  orbits  and  they  are  also  initially  on  the  frame's 
negative  X-axis.  The  sun  is  placed  at  one  A.  U.  with  a 
circular  velocity  at  that  point  solely  in  the  negative 
Y-direction.  The  sun's  X-component  is  transformed  into  the 
correct  units  and  the  same  with  the  velocity.  The  moon's 
position  is  also  on  the  negative  X-axis  at  one  moon  distance 
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The  moon's  position  and  velocity  vector  are  also  given  by: 


• 

<“  i . 

0. 

r  « 

0. 

r  = 

-.229971 

m 

m 

0. 

0. 

EQ  11-16 


The  colony's  initial  position  and  velocity  vectors  are 
given  by  Ref  [7:62].  However,  Ref  [3:31]  provides  the 
vectors  in  an  earth-centered  frame.  The  correction  that 
Beekman  made  is  sufficient  for  the  position  vector,  but  the 
velocity  vector  has  the  wrong  units.  The  units  used  are  in 
MD/TU,  mean  moon  distance  per  absolute  time  unit,  while  the 
units  required  are  in  MD/DAY,  mean  moon  distance  per  day. 
Ref  [3:22]  provides  the  TU  relationship  and  this  is  used  in 
the  correction.  Finally  the  position  and  velocity  vectors 
are  given  by: 
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The  effect  of  the  remainder  of  the  solar  system  on 

the  satellite's  motion  needs  to  he  investigated.  Since  the 
motion  of  the  sun  and  moon  are  known,  the  only  terms  that 
need  to  be  considered  in  the  "n-body"  equation  are  the  motion 
of  the  Earth  and  that  of  the  satellite.  Even  though  this  is 
an  "n-body"  problem  most  of  the  mass  of  the  solar  system  is 
at  such  a  distance  that  the  effects  are  going  to  be  minimal. 
With  that  in  mind,  the  problem  can  be  reduced  to  the  effect 
of  Jupiter's  tidal  acceleration  on  the  Earth's  and 
satellite's  motion.  Jupiter  is  the  closest  of  the  massive 
planets  and  would  reflect  the  largest  effect  on  their  motion. 
As  a  comparison,  the  tidal  acceleration  of  Jupiter  can  be 
related  to  me  sun's  tidal  acceleration.  The  tidal 
acceleration  equation  can  be  written  as: 


“cJ-Gmj(%jreJ 


-  -r  r“3) 
cJ  c  J 


EQ  11-18 


where  rcj“rej+rc*  Expanding  the  cubic  term  of  the  colony: 
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where  r  ?  c  a  a  b  e  r  ix  • .  s  r  d  ^  ! 

<■  i* 
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ininc  terms  n 


-  3  -  3  -  3 

r  =  r  -  3 r  (r  •  r  )  +()  EQ  11-20 

c  J  e  J  •  •  J  e  J  c 


Substituting  this  back  into  EQ  IX-18  to  form: 


u  - 

r  -  -Gm  r  r 

cJ.;  J  eJ  eJ 


( 


r  +r  )  (  r 
e  J  c 


-3 

eJ 


-5 

3r  ,(?  i  * 
e  J  e  J 


* 

)) 

c 


Collecting  terms  finally  gives: 


-Gin 


3(r  .  r  ) ( r  +r  ) 
e  J  c  e  J  c 


-  r 


e  J 


EQ  11-21 


Ref  [2:360]  provides  a  relative  value  of  Jupiter's  mass  and 
distance  from  the  selected  frame*  These  can  be  converted 
into  values  relative  to  the  sun  for  the  comparison. 
Replacing  with  its  equivalent,  .001ms,  and  rgj  with 

4.20r  ,  into  EQ  11-21  will  allow  a  numerical  evaluation  of 

s 

the  acceleration.  Noting  that  rej+rc  “  rej  “4.20  rg  ,  EQ 
11-21  can  be  reduced  to: 
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which  demonstrates  the  tidal  acceleration  due  to  Jupiter  or 

the  remainder  of  the  solar  system  is  approximately  10  7  less 
than  that  of  the  sun's  tidal  acceleration.  Therefore,  we 
need  not  concern  ourselves  with  any  solar  system 
perturbations  other  than  the  sun's. 

The  Four  Body  Problem 

The  initial  conditions  of  the  four  body  problem  are 
the  basis  on  which  the  whole  stability  question  lies.  If  the 
wrong  initial  conditions  are  chosen  the  search  for  stability 
will  be  long  and  tedious.  Capt  Wheeler's  model  however, 
demonstrates  linear  stability  in  a  is  frame.  Therefore  the 
initial  conditions  selected  should  have  a  greater  chance  of 
three  dimensional  stability  than  any  selected  by  random. 
Using  the  position  and  velocity  vector  determined  in  the 
analysis  of  the  Wheeler  model  will  provide  the  needed  vectors 
in  two  dimensions.  Since  the  velocity  and  position  vectors 
are  only  in  two  dimensions  we  need  to  add  the  third  dimension 
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these  starting  conditions.  There  tore  it  will  oecorce 

necessary  in  the  a.u  lysis  to  atoait  y  t  ii  e  s  e  conditions  in  tne 

search  for  stability.  The  manner  in  which  the  modification 
is  made  does  not  matter  and  therefore  can  be  done  by  any 
method  available.  Once  the  modification  is  made  the 
equations  pf  motion  can  then  be  integrated  forward  for  a 
relatively  short  period  of  time.  If  the  colony  is  still  in 
the  vicinity  of  L4,  they  can  be  integrated  forward  for  a 
greater  length  of  time.  This  can  be  repeated  until  the 
period  of  stability  reaches  that  as  desired.  Using  short 
integration  time  spans  in  the  initial  search  for  stability 
will  greatly  reduce  the  amount  of  execution  ti  required  far 
the  integration  of  the  eighteen  equations  of  notion. 
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The  span  of  the  ephcneris  was  selected  fa  be 
approximately  one  year.  This  particular  length  was  chosen 
due  to  its  amount  of  integration  necessary  for  the 

generation.  Longer  time  spans  would  be  nice  but  the 
execution  time  increases  linearly  for  the  integrations 
required.  The  span  does  allow  a  sufficient  number  of 
starting  points,  new  moons,  for  the  analysis.  The  addition 
of  the  dynamic  differential  corrector  for  the  velocity 
vectors  of  the  moon  and  sun  produced  a  positional  error 
related  to  time  of  approximately  six  hours  for  both  the  moon 
and  sun  for  the  slightly  less  than  one  year  generation  span. 
This  error  is  slightly  less  than  .07  percent,  or  less  than 
the  error  required  at  the  outset  of  the  problem.  If  the 
error  were  cumulative  throughout  the  study,  then  the  total 
l:  r  r  o  r  f  c  r  a  fit'--  ,  i  .  r  i  .1  wo  u  •  1  a  p  p  r  ;  i .  r  p  w  r  •  :  •_  ,  or 

12.5  days.  This  total  error  would  have  little,  if  any, 
effect  on  the  question  of  stability. 

The  sun's  velocity  was  corrected  three  times,  the 
first  time  for  a  period  of  one  day,  then  ten  days,  twenty 
days,  and  finally  fifty  days.  The  moon  on  the  other  hand. 
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the  titty  day  ret  ere  11  co  period  proa  need  cue  accuracy  required 

and  such  an  investigation  is  beyond  the  scope  of  this  report. 
The  error  in  the  resultant  ephemeris  is  linear  throughout  the 
entire  span.  The  major  cause  for  this  would  be  a  minor  error 
still  exlstant  in  the  initial  conditions  used  in  its 
generation.  The  use  of  a  tabulated  ephemeris  that  contained 
the  entire  span  needed  would  be  of  some  help,  provided  that 
the  accuracy  is  greater  than  that  achieved  by  this  method. 
However,  as  this  analysis  indicates,  such  an  ephemeris  is  not 
needed.  Table  I  provides  selected  portions  of  the  moon  and 
sun  ephemeris  in  ephemeris  coordinates  which  can  easily  be 
compared  to  Ref  8  to  note  their  accuracy.  Table  II  is  the 
state  vector  of  the  sun-moon  conditions  on  l.u  jail  1979. 


The  Wheeler  Orbit 

The  initial  conditions  derived  from  Ref  [3:31]  were 
initially  integrated  forward  for  a  period  of  one  lunar 
synodic  month.  The  data  obtained  in  the  nonrotating  frame 
was  transformed  into  a  rotating  frame  in  which  the  moon 
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Fig  5.  Restricted  Wheeler  Orbit — One  Month 
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Fig  7.  Restricted  Wheeler  Orbit— Twelve  Months 
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Four-Body  Analys is 

The  selection  of  the  Initial  conditions  were  highly 
Influenced  by  the  periodic  orbit  developed  by  Capt  Wheeler  in 
Ref  7.  If  a  stable  orbit  does  exist,  then  the  initial 
conditions  of  such  an  orbit  should  be  close  to  those  given  by 
Ref  [7:62].  The  third  dimension  of  the  vector  can  be  related 
to  the  moon's  vector  by  giving  the  satellite  the  same  z  and  z 
magnitudes  as  the  moon.  The  moon  and  sun's  initial 
conditions  are  obtained  from  ephemeris  generation  techniques 
discussed  earlier.  However,  Ref  [7:62]  requires  the  moon  and 
sun  bath  be  on  th<-  no,:,  .it  i  a  x-u:.  l  s  tor  those  par  Lie  alar 
initial  conditions.  This  condition  occurs  once  every  new 
moon.  Ref  [8:3]  provides  a  list  of  dates  of  new  moon 
occurences  accurate  to  the  nearest  minute.  The  first 
occurance  after  5  Jan  1979  is  that  of  28.263888  Jan  1979. 
The  equations  of  motion  are  then  integrated  to  this  date, 
neglecting  the  satellite's  motion.  Once  the  integration  is 
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the  order  of  4.5  meters  per  second  or  .001  MD/DAY  should  have 
the  desired  effect. 

A  main  program  was  written  which  allowed  the  user  to 
add  these  delta  velocities  and  specify  the  length  of  the 
integration.  The  integration  period  should  be  variable  in 
length  to  aid  in  the  search  for  stability.  The  length  of 
execution  time  of  the  integration  requires  short  integration 
spans  until  a  likely  candidate  for  stability  is  found.  The 
program  was  written  for  an  interactive  user  which  allowed 
quick  observation  and  interpretation  of  results.  Appendix  E 
contains  the  main  program  used  in  this  effort.  Appendix  F 

C* 

contains  the  subroutine  used  in  the  addition  of  the  delta 
velocities  to  the  state  vector.  First,  the  state  vector 

elements  of  the  satellite  are  rotated  to  the  rotating  frame. 
Then  the  delta  velocities  are  added  to  the  rotated  elements 
and  these  are  then  transformed  back  into  the  nonrotating 
frame.  Once  the  initial  conditions  are  those  which  are 
desired,  the  integration  is  executed  for  the  desired  time 
span.  The  time  span  is  split  into  time  steps  which  are 
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.satellite  elements  used  in  tile  'diieelcr  medel  .we  contained  i  r. 
labie  III.  i'ig's  S ,  4,  10,  ana  i  l  are  trajectories  obtained 
for  periods  of  one,  three,  twelve,  and  sixty  synodic  months. 
Adding  a  delta  velocity  of  .001  M D/D AY  to  all  of  the 
satellite's  velocity  components  of  Table  III,  produces  the 
state  vector  contained  in  Table  IV.  The  trajectory  plots  for 
one,  three,  twelve,  and  sixty  synodic  months  are  contained  in 
Fig's  12  through  15.  Adding  a  delta  velocity  of  -.001  MD/DAY 
to  each  of  the  velocity  components  in  Table  III,  produces  the 
state  vector  contained  in  Table  V.  The  resultant  trajectory 
is  plotted  for  the  same  periods  as  before  and  are  shown  in 
Fig's  16  through  19. 

The  question  of  stability  is  of  paramount  importance 
to  the  problem.  Linear  or  periodic  stability  does  not  exist 
in  the  three  dimensional  model,  but  stability  must  be 
determined  and  in  an  easily  observable  fashion.  To  that 
effect,  a  cross  section  of  the  orbit  is  obtained  by  slicing 
through  the  Lagrangian  point,  using  the  yz-plane.  For  long 
term  stability  to  occur,  the  cross  sections  should  fill  in 
separate,  but  definite,  areas  on  the  plane.  The  reasoning 
behind  this  is  if  an  orbit  is  stable  then  after  a  reasonable 


32 


Tab  K>  Lit.  Slate  Vm:  tor  a  i  J  S  .  2  h  i ;5  !•»  S  .)  ,j  n  1  '■)  j 
i  it  c  lot!  l  n  ^  ..no  i  I  t  i"  . .  1  .,i  at  s 

I  Ml  Tit'll  Cl  It1 1 1  I  i  1 1  it)'  i  Cl  C  I  ilk 

•:  (I  •  : .  :  I  I '»  vr. : 

A  i  ,i  -  ,»o  ! .  •'  1  /  .  ' 

x  1 3  i  •• ,  o  i  •: !■; 7  i o 7 ;< 5 a  /  s? 

X  <  4  )  -  5  .433476 266665 
y  <:;  >  a  .  oh  t/./mh:1:1  i  a 
X  <  6  )  -  .  000 1 928591 60 1 II 4 5 
X  <  7  >  =  .  5701945487 2 2 5 
X<8>  =  -.7312694719864 
X  <  9  >  -  .05406587955571 
X<10)  =  .1940024128416 
X<11>=  .1502478807353 
X< 12)=  -.01690467092801 
X<13)=  -.1904858617353 
X< 141=  -1.032237274618 
X< 15 ) -  .05406587955571 
X  <  1 6 )  -  .  2008356 1 32464 
X ( 1 7 ) -  -.059826122177 
X  < 18 ) =  -.01690467092801 

ROTATING  FRAME  COLONY  POSITION  AND  VELOCITY 
X=  -.7363274298999  XD0T=  -.1706730983 
Y*  .8156863968899  YD0T=  -.121592771 
Z*  .05406587955571  ZH0T=  -.01690467092801 
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Pig  9»  Table  III  Wheeler  Orbit — Three  Months 
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Table  IV 


State  Vector  ot  1 S  .  J  i  >  S I  S  S  Jan  1  J  /  ') 


with  Wheeler  K  I  ■■wont  ••  i  ■  a.' 


nil.  i  i  iii  i 'i  um  1. 1  i  n n-,  Mi  i :  i  niv 

X  (  1  )  23'’  *  |  'V, 0  I  1  .  ' * 1  !'■ 

)  ■■  .  i  ■  '■ ' 

;■(  (  3 )  ■■■■■  .  O  I 

X  (  4  )  -  5.43347626 . V. . 

X  (  5  )  -  4 . OO35534022  I  4 

x <  A  > —  .000 1  92059  i  M'<  i  ;:-r. . 

X<7>  =  .5701 945407:*  I'!; 

X<8)  =  - . 7312694719864 
X<9)~  .05406587955571 
X<10>=  .1940024120416 
X<11>=  .1502478807353 
X  <  12 )  =  -.01690467092001 
X  < 13 )  =  -.1904858617353 
X<14>=  -1.082237274618 
X  ( 15)—  .05406507955571 
X  < 1 A  )  -  .1994321000  19  1 

X< 17 ) -  •  .05965241010420 
X  < 1 8 ) =  - . 0159046709200 1 

ROTATING  FRAME  COLONY  POSITION  ANH  VELOCITY 
X—  -.7363274290999  X)iOT=  -.1696730903 

Y=  .8156863968099  YD0T=  -.120592771 
Z-  .05406587955571  ZD0T=  -.01590467092801 
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Fig  15.  Table  IV  Modified  Wheeler  Orbit  —  Sixty  Months 
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Tab  1 


Still'  '.''I'tiir  dI  S  S  .  J  (>  j  -  I  i  n  !'•/')  wit 

!.v ! i  ,•  ■  >  i.  .•  r  il  I  iii  -  it  i  •>  :  i  ;  .i  ,  ,  „• 


I  i  J  I  I  I  in  i  i  ii  'll  |  I  I  i  i,  i-  .  :  II  lie 


'(<  1  >  ■ 

7  1  1  n  |  .  • 

'  1.  "  !  ■■ 

<  (  ; ;  ) 

,  ( >  j  i! »  o  ,  •:  ,* 

X  (  4  )  - 

X  (  5  ) 

4.00355348221  4 

X  <  6  )  - 

.  OOOI  77059160  1  >145 

X  (  7  )  - 

.5701945407225 

X  <  8  )  = 

-.7312694719064 

X  <  9  )  = 

.05406587955571 

X(10)  = 

.1940024128416 

X(ll)  = 

.1502478807353 

X<12)= 

-.01690467092801 

X< 13 )= 

-.1904858617353 

X<  14 )  = 

-1.082237274618 

X<15)- 

.05406507955571 

X(  L  A  )  - 

.2077391 104736 

X  <  1 7  )  = 

-.05999902616972 

X  <  18  >  = 

-.0179046/092001 

ROTATING  FRAME  COLONY  POSITION  ANIi  VELOCITY 
X—  -.7363274298999  XDOT-  -.1716730703 
Y=  .0156063768099  YDOT^  -.122592771 
Z=  .05406507955571  ZHOT=  -.01790467092001 
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Fig  17.  Table  V  Modified  Wheeler  Orbit — Three  Months 
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i  n  r  i  i;  10.  11  r. 

i  i  i  a  r  i 

/  ,  c  h  e  orbit  s  -  r  o 

rig's  l;  ar.u  '  1  '  • 

a  r  v 

analyzed  and 

their 

cross  sections  are 

presented  in  r'  i  g  ' 

s  21 

ana  12  respect! 

v  e  1  y  • 

An  analysis  of 

tlie  plots  show 

the 

orbits  to  be  unstable  for  all  three  cases. 

The  search  was  continued  in  the  manner  described  above 
until  a  likely  candidate  was  found.  The  plots  for  the 
system  are  contained  in  Fig's  23  through  30.  Fig  30  is  the 
plot  of  the  cross  section  for  a  period  of  fifty  years.  The 
cross  sections  fill  only  the  definite  areas  inscribed  and  do 
not  fall  out  these  areas.  This  leads  to  the  conclusion  that 
this  particular  state  vector  produces  a  stable  trajectory  for 
the  require!  tire. 

The  state  vector  contained  in  Table  VI  is  not 
coincident  with  a  new  moon.  Indeed,  it  is  5.01  days  prior  to 
the  new  moon,  but  an  important  point  is  the  initial  position 
is  precisely  the  initial  position  of  the  Wheeler  orbit.  The 
velocities  differ  only  in  the  third  decimal  place.  The  state 
vector  is  integrated  forward  to  the  new  moon  to  find  the  set 
of  Initial  conditions  required  for  a  stable  orbit.  This 
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Fig  20.  Cross  Section  of  Orbit  of  Fig  11 
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■  :  '  j  .  1  t  !  '>  1 ')  !  n  c  !  -i 

•  •  •  ’  ■  i 

I  jvi  f  r  I  I  T  )  •  ( tJ 1  ‘  "M  '  I  1 1 1 •• 

' 1  I  i  :  1 

,U2>-  V'  :  .  •)  !  V  • 

X  (  X )  ■■■■■■  » o  I '  .48  I  4::44/.9:-U 

X  (  4  >  =  5  .  7"?':>97,A<?004X4 
X  (  S  >  '■  .  !  .94  /VI!! /V  I  ! ! 

X ( A ) =  . 0001 6949789925  I  5 
X<7)=  -,5624053120187 
X  <8 ) =  -.7970584471 
X<9>  =  .08042470202791 
X(10>=  .1985222851766 
X<11>*  -.120556823555 
X< 12)=  .007137337355667 
X< 13)=  -1.090992584716 
X  < 1 4 ) =  -.1313680409769 
X  <  I  5  )  -  .  O  7 ; !  4  7  3  I  5  * ,  I  1 7  0 i ! 

X ( 1 6 ) —  . 01046712 1059 18 
X  < 1 7  >  =  -.2064773418006 
X  < 1 8 ) =  .01 8262800 1 0529 

ROTATING  I" KAMI!  f.Ol ONY  POSITION  AMU  Oi  l  OCT  I 
X=  -.7363274299  XD07--  -.162673098 
Y*  .81568639689  YDOT*  -.127592776 
Z»  .07847315514708  ZliOT=>  .01826280010529 
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Fig  24.  Table  VI  Stable  Orbit  Candidate--Three  Months 
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Fig  26.  Table  VI  Stable  Orbit  Candidate  —  Sixty  Months 
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Fig  27.  Cross  Section  of  Orbit  of  Fig  26 


57 


.  60 


OQ  *2 

09*1 

QZ'l  OQ’O 

0*’Q 

0Q*0‘ 

SIXU-A 

YX-flXIS 


VIII  contain.-.  the  .state  vector  o:  13.  I1  .Ian  1  1  }  w  .l.i  toe- 
stab  I  e  orbit  elements. The  subsequent  plots  and  cross  sections 
are  contained  in  rig's  31  through  j  v •  Inese  plots  are 

essentially  the  same  as  the  plots  that  are  produced  from 
Table  VI,  and  the  conclusion  can  be  reached  that  the  state 
vector  is  at  least  marginally  stable.  Further  investigation 
of  this  state  vector  does  indeed  show  marginal  stability  for 
the  system,  but  the  results  are  inseparable  from  those 
produced  from  Table  VI,  and  are  therefore  not  reproduced 
here . 

Finally,  Ref  [3:54]  refers  to  a  180°  out  of  phase 
orbit:  as  proposed  Fe  f  £  .  To  so.:  re-  f  -  r  this  i:  the 

state  vector  iron  Table  VI  was  integrated  tor ward  tor  tit  teen 

time  steps,  or  half  the  period.  The  satellite's  position  and 
velocity  vector  is  then  recorded  for  use  in  the  system 
analysis.  The  state  vector  of  Table  VI  is  then  combined  with 
this  position  and  velocity  vector  and  tested  for  stability. 
The  vector  produces  a  similiar  orbit  to  that  produced  in 
Table  VI  and  is  reproduced  in  Table  IX.  The  orbits  and  their 
cross  section  for  a  period  of  sixty  months  are  produced  in 
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Tibi.'  VI  r.  still,-  Vector  of  2S.  2(»3°88  .Inn  1979  [n<- lu.!  i  n:- 

S  t  1 1>  I  —  foil  1  ■  .  ■  . '  ■ 

i:  f  J  I  T  I  r.t  i  i:  111  I  i  iO.lli  1 1  I  UP 
X  ( .1.  )•••'  ■  :  i  7/>o  i  :  . 

/.  ( ?  ,<>  :  .  '  :  '  ' 

X  ( 3 ) -  .  ( ■  I  ’ ; 7 3 0 7 "V /; 7  '■ 

X  ( 4  !  •-  !.)  I-  1 . > , > 4  ■  •  >.' (i 6 <!> 6 1 1 
X  (  5  >  -  4.08355 34133"'!  4 
X  (  6  >  •»  .  0 00  I  7 77'  . 7 .1  X 0 1  ft 9 •/, 

X(7  >  =  .5701945487225 
X  <8 )  =  - .  7312694719064 
X  <  9 ) =  .05406587955571 
X(10>=  .1940024128416 
X(ll)=  .1502478807353 
X<  12>=*  -.01690467092801 
X< 13 )=  -.5646495963434 
X<14)=  -.9507939452589 
X<15)=  .1214735632235 
X  < 1 6 ) =  .181 5202065839 
X< 17 ) ~  -.09474670 2 2  5  I  7  4 
X ( J. 8 >  =  -.002393534008842 

ROTATING  FRAME  COLONY  POSITION  AND  OFL..OCITY 
X»  -.4025971030472  XDOT  =  -.1863345732500 
Y®  1.029929058532  YDOT®  -.08488786601906 
2“  .1214735632235  ZDOT=  -.002393534000842 


TT 
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Table  Via.  Sons  i  L  i  /  i  t  y  Tost  State  Vector 

INI  l  I  .  .I  CONIi  l  i  I  Dili;  V'l  i.  I  UK 

;  j  ■  : ) !  ■  -; V  •*,</  , 

,'3o'. 

X  ( 3 )  ~  -  .  o  1 55242696305 
X  <  4  )  -  5 . 700942315224 
X  (  5 )  «  3 . 5  /,  ?  9 09  4  7 0  A  9 
X  <  6  >  ~  . 0001 602791 094 1 72 
X<7>  =  - « 6116523585634 
X<8 )  =  -.7651139719590 
X<9)  =  .07847315515033 
X<10)=  .1901730539485 
X<11)«  -.1315409385366 
X(12)«  .008262800093982 
X< 13 ) =  -1.096901223061 
X< 14 ) ®  -.06580341132929 
X  < IS  1  =  . 07847315514700 
X< 16)^  -.001915574420106 

X  <  1  7  >  -  - . 2047336060609 
X< 18 ) *  .01826280010579 

ROTATING  FKAME  COLONY  POSITION  AND  VELOCITY 
X*  -.7363274299  XIHIT-  -.162673090 

Y-  .81568639689  YDOT-  -.127592776 
2»  .07847315514708  ZD0T»  .01826280010529 
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rig  32.  Table  VIII  Sensitivity  Orbit — Three  Months 
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Pig  35.  Cross  Section  of  Orbit  of  Fig  34 
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Pig  37.  Cross  Section  of  Orbit  of  Fig  36 
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Fig  39.  Cross  Section  of  Orbit  of  Fig  38 
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Table  IX 


State  Vector  of  180°  Out  of  Phase  Orbit 

INITIAL  CONDITIONS  VECTOR 
X<  1  )  =  204.3194374608 
X(2)=  -324.0192778463 
X<3)  =  -.01548148447034 
X(4)  =  5.772936900824 
X<5)  =  3.594798567616 
X  ( 6 )  =  .000169497899237 
X ( 7 ) =  -.5624053126209 
X  <  8 )  =  -.7970584467333 
X  ( 9 )  =  .08042470200617 
X<10)=  .198522285084 
X(ll)=  -.120556823687 
X(12>=  .007137337369057 
X  <  1 3 )  =  -.8200582334128 
X<14)=  .4868110858079 
X< 15)-  ,0785 
X<16)=  -.07662306948368 
X< 17)=  -.2249784474187 
X< 18)=  .0213 

ROTATING  FRAME  COLONY  POSITION  AND  VELOCITY 
X=  -.0750263  XDOT=  -.228 
Y=  .950711099  YDOT=  -.0670999 
Z=  .0785  ZDOT  =  ,0213 
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Fig  40.  180°  Out  of  Phase  Orbit--One  Month 
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Fig  43.  180°  Out  of  Phase  Orbit--Sixty  Months 
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Fig  44.  Cross  Section  of  Orbit  of  Fig  43 
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Fig  48.  Cross  Section  of  Orbit  of  Fig  47 
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Pig  51.  Cross  Section  of  Orbit  of  Fig  50 
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The  primary  conclusion  of  the  report  is  stable  orbits 
do  exist  in  the  three  dimensional  system.  These  orbits  are 
stable  for  at  least  fifty  years  and  allow  planners  to  select 
vehicles  whose  lifetimes  compare  with  this  period  of 
stability.  Additionally,  The  stable  orbits  are  nominally 
insensitive  to  injection  error  and  can  be  maintained  with  a 
small  initial  correction  to  the  actual  stable  orbit.  The 
addition  of  a  marginally  stable  orbit  180°  out  of  phase  with 
the  stable  orbit  increases  the  coverage  available  to  remote 
sensing  equipment  stationed  in  the  L.irraiv  i.in  vicinities. 

twenty  years,  which  is  a  sufficient  period  for  an  unmanned 
satellite.  The  use  of  a  controller  would  probably  insure  the 
stability  of  that  particular  orbit  for  an  additional  ten 
years  . 
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The  stable  orbit 

wo  u  1  d 

not  have 

been 

found  if  the  error  had  not  occurred.  Further  ,  from  an 
analysis  of  the  pertinent  figures,  if  a  different  ephemeris 
was  selected,  the  orbit  would  not  have  been  found  either* 
Wheeler's  initial  conditions  are  only  approximately  matched 
in  the  first  three  months,  and  then  again  only  after  a  period 
in  excess  of  five  years. 

The  error  and  ephemeris  selection  were  tested  to 
define  the  uniqueness  of  the  stable  orbit*  The  state  vector 
of  Table  VI  was  integrated  forward  for  a  random  length  of 
time,  where  the  '.'heeler  conditions  were  inserted  in  the 
prospect  of  reproducing  the  stable  orbit.  While  not 
conclusive,  the  results  indicated  that  the  conditions  were 
unique  for  a  short  time  period,  i.e.  five  years.  The  fact 
that  a  stable  orbit  using  the  Wheeler  conditions  at  a  new 
moon  was  not  discovered,  does  not  discount  the  possibility  of 
one  existing.  The  fact  that  the  orbits  are  very  ephemeris 
oriented  should  allow  for  the  discovery  of  one  using  the 
correct  ephemeris.  The  search  for  such  an  orbit  was  not 
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of  controllers  on  all  orbits  described,  to  increase  their 
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orbit  using  the  Wheeler  initial  conditions,  should  be 
continued  to  prove  or  disprove  the  possibility  of  a  stable 
orbit  existing  with  a  starting  point  as  a  new  noon.  Finally, 
the  orbits  described  should  be  extended  to  the  end  of  their 
stability  in  a  search  of  the  stable  lifetime.  The  orbit  in 
Table  VI  was  actually  integrated  in  excess  of  sixty  years  and 
showed  no  sign  of  decaying  into  instability. 
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Appendix  A 

Subroutine  Containing  the  Equations  of  Motion 


SUBROUTINE  F  <  T»  X  *  DX» NEQN ) 

COMhON/DFIL/GSUN r  MBA  » MDE *  MS  t MU » GE  >  PI 
DIMENSION  X(18)rDX(18) 

REAL  MDE* MS  >MDA* MU 

RS=SQRT( <ABS(X(1) ) ) **2+ < ABS ( X < 2 ) ) >**2+(ABS(X<3) > >**2> 
RM=SQRT ( <ABS<X(7) ) )  **2+ < ABS < X < 8 ) ) > **2 M ABS ( X < 9 > ) )*#2) 
RC=SQRT( <ABS(X(13) ) ) **24 ( ABS ( X ( 14 ) ) ) **2+ ( ABS< X ( 15 ) ) >**2> 
RCS==SQRT  ( <  ABS<  X  ( 13 ) -X  ( 1 )  )  )*#2+< ABS<X< 14)-X(2) )  ) **2+ ( ABS<  X ( 1 
&5)-X<3) ) )**2) 

RCM-SQRT (  (ABS(X(  13)-X<7)  )  ) *#2  +  ( ABS (  X ( 14 ) -X ( 8 ) ) ) **2+  ( ABS ( X ( 1 
&5>-X<9> ) >**2> 

RMS=SGRT < (ABS(X<7)-X< 1 > ) ) **2+ ( ABS ( X < 8 > -X ( 2> >  >**2+<ABS(X(9 >- 
&X<3) ) )**2) 

DX ( 1 ) =X ( 4  > 

DX<2)=X<5) 

DX<3)=X(6) 

DX  <  4  )  ---GSUN>KX  ( 1 )  /RS**3 
DX ( 5 ) =-GSUN*X ( 2 ) /RS**3 
DX ( 6 ) =~GSUN*X ( 3 ) /RS**3 
DX ( 7) =X ( 10 ) 

DX(8)=X  < 1 1 ) 

DX(9>=X( 12) 

DX< 10)=GE*(-X(7)/RM**3-MS*<  <X<7)-X( 1) )/RMS**3+X( 1 >/RS**3> ) 
DX<11)=GE*<-X<8)/RM**3-MS*<  <X<8)-X<2) )/RMS**3+X(2)/RS**3>  > 
DX<12)-GE*<~X<9)/RM**3-MS*<X(9)-X<3) > /RM3**3+X ( 3 ) /RS**3 ) 

DX< 13)=X( 16) 

DX( 14 ) =X  < 17 ) 

DX ( 15 ) =X ( 18 ) 

DX(  1A)=GE*(-<1-MU)*X<  13)/RC*#3-MS*<  <X(  13)-X<  1 ) >/RCS**3 
S+X<1)/RS**3)-MU*< (X< 13)-X<7) )/RCM**3fX<7)/RM**3) ) 
DX<17)=GE*<-<1-MU)*X< 14>/RC**3-MS*<  ( X ( 14  ) -X  ( 2  )  )  /RCS**3 
&+X(2)/RS**3)-MU*< (X< 14)-X<8) ) /RCM#*3+X ( 8 ) /RM**3 ) ) 

DX<  18)-~GE*<-<  1-MU)*X<  15>/RC**3-MS*(  (X(  15>-X<3>  >/RCS**3 
&+X(3)/RS**3)-MlJ*<  <X(  15)-X<9)  )/RCM**3+X<9)/RM*#3>  ) 

RETURN 

END 
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Appendix  B 


Subroutine  For  Determining  the  Velocity  Correction  Vector 
SUBROUTINE  ENT ( F * TOUT * I 1 , 12 * A1 * A2 * A3 , X ) 

COMMON  /DFIL/GSUN * MBA » MDE * MS * MU » GE » PI 
DIMENSION  X<18) »XC<18> *XD<18> 

DIMENSION  R1 (3) *AC<3) 

DIMENSION  A(3*3)»B<3>»C(3*3) 

DIMENSION  WKAREA <800) * IW0RK<30) 

REAL  MS*MU*MDE*MDA 
REAL  LAT  *LON 
DO  10  NK=1 * 3 
DO  1  1=1*18 
XC<I)=X<I) 

XD< I >=0. 

1  CONTINUE 
IF(I1.GT,1)G0T09 
AC< 1 )=A1 
AC<2)=A2 
AC<3)=A3 

9  IF<I1.LT.7)G0T011 

AC<1)=A1*COS<A2)*COS<A3> 

AC<2)=A1*SIN<A2)*C0S<A3) 

AC<3)=A1*SIN<A3) 

11  NEGN=18 

ERR-1E-9 
IFl.AG=l 
T=0 « 

CALL  ODE ( F  *  NEON rXC,Tr TOUT  *  ERR  *  ERR *  I FLAG  *  WKAREA *  I WORK ) 
DO  2  1=1*18 
XD ( I ) =XC ( I ) 

2  CONTINUE 
13=11+2 

DO  3  1=11*13 

R1(I- (11-1))=  <XC( I >-AC( I-(Il-l))) 

3  B(I-(I1-1) ) =XC ( I ) 

DELTAV=lE-4 

DO  4  1=1*3 
DO  8  K=1 *  18 
8  XC(K)=X(K) 

XC < 1+ < 12-1 > )=XC<I+ <12-1 ) )+DELTAV 
IFLAG=1 
T=0. 

CALL  ODE(F*NEQN*XC*T »TOUT * ERR* ERR* IFLAG* WKAREA* I WORK) 
DO  5  J=1 *  3 

A< J»I)=XC< J+<I1-1) )-XD< J+<I1-1)  ) 

A<J»I)=A<J*I) /DELTAV 
CONTINUE 
CONTINUE 
IDGT=7 
N=3 
Nl  =  l 

CALL  LIN02F <A*N*N*C*IDGT * WKAREA* IER ) 

CALL  VMULFF <C*R1*N*N*N1*N*N*B*N*IER) 

DO  7  1=1*3 

X<I+(I2-1 ) )=X< I+<I2-1 ) )-B< I) 


91 


Appendix  C 


Subroutine  Used  to  Transform  to  the  Rotating  Frame 


SUBROUTINE  COLROT ( I , N , T , X , PI > 

DIMENSION  X ( 18 ) 

THETA=ATAU2(X<8>  ,X(7) )+PI 

IF ( THETA . LT . 0 ) THETA=THETA+2 , *PI 

A=X (13) *COS ( THETA ) +X  < 1 4 ) *S I N ( THETA ) 

B=-X ( 1 3 ) *S IN ( THETA ) +X ( 1 4 ) *COS ( THETA ) 

C=X (15) 

AD=X( 16 ) *COS( THETA )+X( 17 )*SIN( THETA) 

BD=-X ( 16 ) *SIN( THETA )+X( 17) *COS( THETA) 

IF(I.NE.N)G0T01 
PRINT* r ‘  * 

PRINT*,*  ’ 

PRINT*f "ROTATING  FRAME  COLONY  POSITION  AND  VELOCITY* 
PRINT*,*  * 

PRINT*, *X=  * , A, *  XDOT=  * , AD 
PRINT*,*  * 

PRINT* , * Y=  * , B , *  YDOT=  *,BD 
PRINT*,*  * 

PRINT*, *Z=  *,C,*  ZDOT=  *,X(18) 

PRINT*,*  * 

WRITE(6)  T, A,B,C 

RETURN 

END 
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Appendix  D 
Plotting  Subroutine 


A=X ( I ) 

B=X ( 1+1 ) 

IF< (A-XL)#<B-XL> ,GE,0)G0  TO  4 
XIN=ABS (X(I)-X(I+1)) 

DX=XL-A 

DX=DX/XIN 

YIN=Y<I+1)-Y(I) 

BY=DX#YIN 

L=L+1 

YX<L)=Y( I )+BY 
ZIN=Z(I+1)-Z(I) 

DZ=DX#ZIN 
ZX<L)=Z(I)+DZ 
4  CONTINUE 

CALL  PLOT (0«»0.r-3) 

PRINT# »  * NUMBER  OF  PLOT?* 

READ  15»HFD 

CALL  SYMBOL ( 0 . 1 0 ♦ t . 25  f HFD , 0 . r  2 ) 

15  FORMAT <1A2) 

CALL  PLOT («5r»5»3> 

CALL  PLOT (9»25t  « 5  f  2 ) 

CALL  PLOT <9.25»A.5»2) 

CALL  PLOT <  «5»6.5»2) 

CALL  PLOT <«5r«5»2) 

CALL  PLOT < 1 . » 1 «  » -3 ) 

CALL  SCALE (YX»7«75»L»1) 

CALL  SCALE<ZX»5. »L» 1 ) 

CALL  AXIS (0 . » 0 . » 7HYX-AXIS»-7 » 7 » 75 f 0 . r YX ( L+l ) , YX < L+2 > ) 
CALL  AXIS <0 . ,0. .7HZX-AXIS»7r5. r90. rZX(L+l ) ,ZX(L+2> ) 
CALL  LINE <YX»ZX»L>l»-2rl) 

CALL  LINE<YX,ZX»L»1>-1»4) 

10  CONTINUE 

PRINT#  f *  HOW  MANY  TO  PLOT?  * 

READ# » N1 

DO  11  11=1 »N1 

XX(II)=X<II) 

YY  < 1 1 )=Y< 1 1 > 

11  CONTINUE 

CALL  PLOT <15.?0«»-3> 

CALL  PLOT <-.5»-«5»3> 

CALL  PLOT <8. 25 5*2) 

CALL  PLOT <S.25»5.5f2) 

CALL  PLOT (-.5r5*5»2) 

CALL  PLOT (-.5»-.5f2) 

CALL  SCALE (XXr7»75»Nlfl> 

CALL  SCALE (YY  »5.»N1»1) 

CALL  AXIS< 0. tO, f6HX-AXIS»-6f7.75fO»  »XX(N1  +  1) » XX  <  N1 +2 ) ) 
CALL  AXIS(0. ,0. rAHY-AXIS»6r5. * 90 . r YY < N1  + 1 >  > YY < N1 f 2 ) ) 
CALL  LINE(XX»YYfNl/lf30.»2.) 
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Appendix  E 


Main  Interactive  Program 


PROGRAM  MA I N ( INPUT » OUTPUT *  TAPES * TAPE6 *  TAPE 7 *  TAPES ) 
DIMENSION  X( 18) * WKAR£A<  800 ) * IU0RK(30) 

COMMON  /DF IL/GSUN » MDA * MDE *  MS  *MU *  GE  *  PI 
REAL  MS * MU  * MDA * MDE 
REAL  LAT*LON 
EXTERNAL  F 

DATA  GSUN*MDA  *MDE *  MS  *  MU  *  GE/1 . 7442438E4  *  2 . 5695187E-3 * 
860.268165*328912.  *  .0121506683* 5.2386349E-2/ 

PI=ACOS ( -1 .  ) 

INITIAL  POSITION  VECTOR 


PRINT# *  * #3  EQUINOX  OF  1950.0" 

REWIND  5 
READ (5)  T*X 
REWIND  6 
NEQN=18 
T=0. 

1  =  1 
N=1 

DEl.TAT=29 . 530589/30 , 

TOUT=DELTAT 

ERR=lE-9 

PRINT** " INPUT  NUMBER  OF  DAYS  TO  INTEGRATE - 1074  OR  LESS" 

READ*»NN 

CALL  COLROT ( I , N * T *X*PI) 

REWIND  6 

CALL  COLLOC ( T  * X * F'I ) 

PRINT**"  " 

CALL  COLROT (I*N*T*X*F'I) 

NN=NN+1 
DO  2  1=2 »NN 
IFLAG=1 

1  CALL  ODE ( F  *  NEQN  *  X  *  T  *  TOUT  *  ERR  *  ERR  * IFLAG  *  WKAREA  * IWORK ) 
T=TOUT 

TOUT=TOUT+DELTAT 

CALL  COLROT <I*NN*T *X*PI) 

2  CONTINUE 
REWIND  7 
WRITE < 7) T  *X 
STOP 

END 
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Appendix  F 

Subroutine  Used  to  Locate  the  Orbit 


SUBROUTINE  COLLOC ( T , X , PI > 

DIMENSION  X< 18 ) 

THETA=ATAN2  (X ( 8) » X( 7 ) )+PI 
IF  (THETA. LT.O)  THETA=THETA+2 . *PI 

PRINT* » ‘INPUT  NEW  INITIAL  CONDITIONS?  1=YES,0=N0' ' 

READ*  ,  N 

IF(N,EQ.1)G0T01 

PRINT*, *  INPUT  DELTAS  TO  INITIAL  CONDITIONS?' 

READ*,  N 

IF(N.EQ.O)  GOT02 

PRINT*, 'INPUT  DELTA  VALUES' 

READ*,  A , B , C , A D , E< D » C D 

A1=X( 13 ) *COS ( THETA )+X( 14) *SIN( THETA) 

Bl=-X ( 13 ) *SIN ( THETA )+X( 14 )*COS( THETA) 

ADI =X( 16 ) *COS ( THETA ) +X( 17) *SIN( THETA) 

BD1=~X( 16 )*SIN( THETA )+X( 17 )*COS( THETA) 

A=A+A1 

B=B+B1 

C=X<15)+C 

AD-AD+AD1 

BD=DD+BD1 

CD=CD+X  < 18 ) 

G0T04 

1  PRINT*, 'INPUT  NEW  INITIAL  CONDITIONS' 

READ*,  A , B , C , AD , BD , CD 

A  CONTINUE 

X ( 13 ) =A*COS< THETA )-B*SIN< THETA) 

X< 14  >=A*SIN (  THETA ) +B*COS (THETA ) 

X ( 15 ) =C 

X ( 16 ) =AD*CGS( THETA )-BD*S IN (THETA) 

X ( 1 7 ) =AD*S IN ( THETA ) +BD*COS ( THETA ) 

X( 18) =CD 

2  PRINT*, 'DO  YOU  WANT  THE  CONDITIONS  ON  PFILE?' 
PRINT*,'  ' 

READ* , N 
REWIND  8 

IF(N.EQ.1)WRITE(8)T,X 
PRINT*,'  ' 

PRINT*,'  INITIAL  CONDITIONS  VECTOR' 

DO  3  1=1,18 
PRINT*, '  ' 

PRINT*,'  X(*,I,')=  * , X ( I ) 

3  CONTINUE 
RETURN 
END 
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